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Analytical nuclear gradients for fully internally contracted complete active space second-order per¬ 
turbation theory (CASPT2) are reported. This implementation has been realized by an automated 
code generator that can handle spin-free formulas for the CASPT2 energy and its derivatives with 
respect to variations of molecular orbitals and reference coefficients. The underlying complete active 
space self-consistent field and the so-called Z-vector equations are solved using density fitting. The 
implementation has been applied to the vertical and adiabatic ionization potentials of the porphin 
molecule to illustrate its capability. 


Substantial effort has been devoted to implementing com¬ 
plex formulas in quantum chemistry, resulting in accurate and 
useful computational tools for chemical applications. How¬ 
ever, some equations are too complicated to be handled man¬ 
ually. Among such examples is the analytical nuclear gradi¬ 
ent theory of fully internally contracted complete active space 
second-order perturbation theories (FIC-CASPT2 or simply 
CASPT2)ii— and its variants," whose complexity has ham¬ 
pered their implementation for more than two decades since 
FIC-CASPT2 was developed.^ In this study we address this 
problem by employing an automatic code generation approach 
to help enable many chemical applications. For instance, such 
implementations can be used for the geometry optimization^ 
of strongly correlated molecules as complex as those investi¬ 
gated by respective single-point calculations. They also have 
the potential to replace mean-field models (e.g., complete ac¬ 
tive space self-consistent field, or CASSCF) often used in 
photodynamics simulations involving the ground and excited 
states of moleculesi^ 

The challenge in implementing the nuclear gradients for 
the FIC-CASPT2 method can be easily recognized as follows. 
The state-specific CASPT2 energy functional for the «-th state 
is 


£ = (d)o|H|Oo) + 2<(I)o|H|<l)i) + (Oil/- £[,"^|Oi>, (1) 

in which / is the standard state-specific Fock operator, and 
- (Oo|/|Oo). The reference and correlated wave func¬ 
tions are defined as (using the reference Cl coefficients c// 
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in which C is the molecular-orbital coefficient matrix param¬ 
eterized by an anti-Hermitian matrix k, 

C = Cinit exp(A:). (5) 


The first term on the right hand side of Eq. (HI is the 
Hellmann-Feynman force^^ and the second term is zero in the 
simplest state-specific CASPT2 case because £ is stationary 
with respect to variation of Tq. The third term also appears 
in the standard single-reference algorithms. The complex¬ 
ity of the equations for FIC-CASPT2 nuclear gradients stems 
from the last term, which is associated with the reference- 
coefficient derivatives; since the first-order wave functions are 
expanded in a basis that is dependent on cf^ [Eq. (O], every 
single term in £ contributes to dSIdcf^ in a nontrivial way. 

The use of partially internally contracted or uncontracted 
basis functionsiS (referred to as WK-CASPT2 in the follow¬ 
ing) for first-order wave functions greatly simplifies the equa¬ 
tions, making it tractable to manually implement nuclear gra¬ 
dients for such variantsJiri^ This is because, in these meth¬ 
ods, parts (or all) of the first-order wave function are expanded 
in terms of excited Slater determinants. 


where I labels Slater determinants, and Eq are external and 
semi-external excitation operators (such as Eai,bj and Ers,at)'^ 
In this article a and b label virtual orbitals, i and j label closed 
orbitals, r, s, t, and u label active orbitals, and x, y, z, and w 
label any orbitals. The energy functional £ is minimized with 
respect to the Tq amplitudes to give the CASPT2 energy. The 
nuclear energy gradients are the total derivative of the energy 
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which are not dependent on c/' (note that in practice only 
distinct determinants are included in the sum). There is also 
an implementation of nuclear gradients of uncontracted mul¬ 
tireference configuration interaction (MRCI)ji^ However, the 





2 


foimal scaling of the size of first-order wave functions in these 
methods is factorial with respect to the number of active or¬ 
bitals, rendering them sub-optimal for large calculations.— 
Over the last decade, various automatic code generation ap¬ 
proaches have been developed to replace tedious, error-prone 
manual implementation processes in quantum chemistryJ^ 
The automated higher-order coupled-cluster (CC) implemen¬ 
tations by Kallay et and by Hiratai2i2£ were the first to 
demonstrate that the automation strategy can produce pro¬ 
grams that are competitive with hand-optimized codes in se¬ 
rial and massively parallel environments, respectively. Han- 
rath and co-workers have realized an arbitrary-order CC code 
generator that is almost optimal— Recently this strategy has 
been extended to various methods such as CC-F12,i22r2^ rela¬ 
tivistic CC^^l local CCr^^ open-shell CCr^ill and excited- 
state CC methods-2i^ Note that these automation techniques 
are complementary to meta-language approaches (such as 
siAL,— iTF,— LiBTENSOR,— TiLEDARRAY,— to name a few), be¬ 
cause code generators can be used to synthesize computer 
codes in any meta language as well. 

Furthermore the automation strategy has been applied to 
development of multireference electron correlation theories. 
Neuscamman et al. used an automated scheme for the canon¬ 
ical transformation theoryj^ Parkhill et al. developed local 
active-space methods (e.g., an active-space perfect quadru¬ 
ples model^) using a sparse frameworkj^ Hanauer and Kohn 
extended their string-based code to implement an internally 
contracted MRCC method-S Saitow et al. reported a fully 
internally contracted MRCI method based on density-matrix- 
renormalization-group reference functions— i 

In this work we extend the automatic code generation ap¬ 
proach to realize FIC-CASPT2 nuclear gradients that have 
been sought for a long time. Following the standard 
approach^iil we use the CASPT2 Lagrangian and the so- 
called Z-vector equation,— instead of directly evaluating 
Eq. (HI, to avoid computation of geometry derivatives of wave- 
function parameters (such as IdR). The state-specific 
CASPT2 Lagrangian isi^ 


£ = £ + itr[Z(A - A')] - ^tr[X(C'^SC -1)] 
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Each term on the right hand side (other than £) corresponds to 
a constraint arising from the CASSCE reference calculation. 
A is the orbital gradient of CASSCE, and Z is its Lagrange 
multiplier. S is the overlap matrix, and X is the Lagrange 
multiplier for orbital orthogonality. The next two terms are 
related to the stationary condition in the full configuration in¬ 
teraction in the active space performed in CASSCE. Here m 
labels states averaged in CASSCE and W,„ is the weight in the 
state averaging. The final term originates from the use of the 
frozen core approximation in CASPT2. See details in Ref.fol 
When the stationary conditions on X. with respect to all the 


parameters and multipliers are met, the nuclear gradients can 
be computed as 

\ 

dR dR /’ ^ 

since only molecular integrals have explicit dependence on 
the nuclear displacement R. We will consider level shifts^ in 
future work, which requires the additional implementation of 
the so-called A equation (as do multistate variants^). 

We have developed a new automated code generator, named 
SMITHS, that derives equations on the basis of the spin-free ver¬ 
sion of Wick’s theorem, in which the normal ordering is de¬ 
fined with respect to the closed-Eock vacuum. In addition, 
SMITHS expresses the terms that involve active-orbital indices 
as a sum of canonical terms so that they can be computed from 
canonical density matrices and its derivatives, e.g., 

(^ol ^ ] t'o-S^^t^pUpl^o} — 2.6rs(J^o\u ~ drt(ro)sit ~ (Fo)^^,?^! 
pa- 

(9a) 

(/| ^ ra-sltlup\(i>Q) = 2b„(Fo)f„ - b«(Fo)j„ - (Fo)^r,r«- (9b) 

pa 

with cr and p labeling spins. Note that (Fo)a = (d>o|£Al®o) 
and (Fq)^ = (/|£a|<1>o) where Ea is a general operator. Here 
SMITHS is used to implement the following expressions: 

(10a) 
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in which |T') = Yji UV) is any multi-configuration reference 
function, and g - YjaGaEa and R. - are general 

and excitation operators, respectively. Note that the determi¬ 
nant index I is treated analogously to the orbital indices in the 
generated code; for instance, (Fo)f„ is viewed as a three-index 
tensor whose size is iidetnact (”det and iiact are the numbers of 
the determinants and the active orbitals, respectively). 

Using this machinery, we automate the nuclear-gradient im¬ 
plementation as follows. First, to optimize Tn the program for 
computing the CASPT2 residual vectors is generated using 
Eq. (llOab . i.e., 

— - 2 [(Q|/- £<"^ 101 ) + (Q|H|d)o)]. (11) 


Second, the Z-CASSCE equation is solved, which is a set of 
coupled equations defined as 
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are implemented by smith3 using Eq. (llOcI) . Next, for 
Eq. ( I12bb the MO-coefficient derivatives of the CASPT2 en¬ 
ergy. 


are calculated from the one- and two-body density matrices 
(see Refs. M and[lj for explicit formulas). The density ma¬ 
trices are defined as 

(ri),^ = 2<(l)o|£,^l'l>i>, (15a) 

(r 2 ).y = - (OolE^^IOoXOilOi), (15b) 

(ri)xy,ZH- = 2 (<Do|£i 3 ,,j„|<t>l). (15c) 

and are implemented using Eq. (llObI) . Given and Yrs, so¬ 
lutions of Z-CASSCE [Z, X and Zij in Eq. (|7]i] can be 
obtained as detailed in Ref.[l^ Using these wave function pa¬ 
rameters and the Lagrange multipliers, effective density ma¬ 
trices are formed, which are then contracted to two-index and 
three-index gradient integral 

The generated code uses a tile-based data layout that is 
similar to those used in xcsi^ and in the earlier version of 
sMiTHi22i2i All the code implemented in the bagel package^ 
and the code generator smith 3^ are openly available under 
the GNU General Public License. The technical details on the 
implementation, working equations, and source code of the 
smith3 program are also found in Supplementary Materials4^ 
CASSCE and Z-CASSCE were manually implemented in 
BAGEL using density fitting (DE) for efficiency as reported in 
Ref.lil In CASPT2, four-index two-external integrals were 
explicitly constructed from DF integrals. The smith3 pro¬ 
gram generated ca. 1150 tasks, the majority of which are 
tensor contractions. Each task is expressed as a node of a 
tree-like directed acyclic graph, which we traverse at runtime. 
This infrastructure should assist in interfacing smith3 code to 
parallel-runtime libraries in the future. 

First, to show the numerical impact of full internal con¬ 
traction on geometrical parameters, we optimized the ground- 
state geometry of the A,A’-diiminato-copper-dioxygen com¬ 
plex [(H 5 C 3 N 2 )Cu 02 ] in its side-on coordination configura¬ 
tion. The ground state is singlet. We used the (14e, 9o) active 
space consisting of an in-plane d orbital of copper and eight 
valence orbitals of dioxygen as suggested in earlier work.-il^ 
The aug-cc-pVDZii^ and def2-QZVPP/JKFITi^ basis sets 
were used for orbital and auxiliary functions, respectively. Ta¬ 
ble |T] compiles optimized Cu-O and 0-0 bond lengths. It is 
evident that neither the degrees of internal contraction (i.e., 
CASPT2 and WK-CASPT2) nor the DF approximation has 
impact on the bond lengths. Our program did not take advan¬ 
tage of spatial symmetry. One optimization step took about 
13 min. using 2 Xeon E5-2650 CPUs (2.0 GHz). More than 
half the time is spent for evaluation of Eq. (fT3]) . Note, how¬ 
ever, that the code generated by smith3 has not been threaded 
efficiently, and there is room for further improvement. 

Next we calculated the vertical and adiabatic ionization po¬ 
tentials (IPs) of the porphin molecule (C 20 H 14 N 4 ) using the 
optimized geometries computed by CASPT2. The cc-pVDZ 


TABLE I. Optimized Cu-O and 0-0 bond lengths (in A) for the 
ground state of (H5C3N2)Cu02 using CASSCF and CASPT2 with 
aug-cc-pVDZ and the (14e, 9<2) active space. DF was used unless 
otherwise stated. 


Method 

Cu-O 

0-0 

CASSCF 

1.886 

1.386 

CASPT2 

1.820 

1.399 

WK-CASPT2_^ 

1.820 

1.400 

WK-CASPT2 (w/o DF)^ 

1.820 

1.400 


“ Partially contracted CASPT2 computed using molpro™ 


basis set^ was used together with the corresponding JKFIT 
basis set^ for DF. The (4e, 4o) and (3e, 4o) active spaces were 
usedr^ which consist of the four frontier orbitals of Gouter- 
man’s model3 The numbers of (correlated) inactive and vir¬ 
tual orbitals were 55 and 323, respectively. One optimization 
step took about 30 min. on the same hardware. For com¬ 
parison, we also computed the IPs from the PBE functional^ 
and MP2 (with an unrestricted variant for the radical cation) 
using TURBOMOLE.— The PBE calculations were performed us¬ 
ing the def2-S VP basis set.— Geometry optimization using all 
methods including CASPT2 was performed without imposing 
spatial symmetry; the geometry of the neutral porphin was 
found to belong to the D 2 /, symmetry group, whereas that of 
the radical cation was found to be € 21 , (even when an initial 
geometry was set to a D 2 ;, structure, optimization converged 
to this minimum). Similar symmetry breaking of metallopor- 
phyrin cation radicals due to the pseudo Jahn-Teller effect 
has been reported in the literature— The optimized geometry 
of the radical cation from unrestricted PBE was found to be 
D 2 / 1 . We were not able to optimize the geometry of the radical 
cation using unrestricted MP2 due to wave function instabil¬ 
ity. The geometrical parameters are compiled in Supplemen¬ 
tary Materials— The IPs are shown in Table lU The vertical 
IP computed by CASPT2 (6.84 eV) is in good agreement with 
the experimental value (6.9 eV)— The difference between the 
vertical and adiabatic IPs computed at 0.18 eV is an order of 
magnitude larger than that computed using PBE. Furthermore 
we computed the IPs with CASPT2 using the PBE-optimized 
geometries. As expected, the vertical IP is almost identical to 
that computed at the CASPT2 geometry; however, the adia¬ 
batic IP is larger than the vertical IP, which attests to the im¬ 
portance of geometry optimization at the CASPT2 level. 

In summary we have used automatic code generation to 
realize analytical CASPT2 nuclear gradients with full inter¬ 
nal contraction. Our implementation has been applied to the 
V,V’-diiminato-copper-dioxygen complex to show that er¬ 
rors due to full internal contraction are marginal. We have 
also computed the vertical and adiabatic IPs of the porphin 
molecule to illustrate the capability of our implementation. 
There is, however, room for improvement in our program in 
terms of efficiency and storage requirement; currently, appli¬ 
cation of our program is limited by the storage of 

(ro)rr',.?5',r/'’ ~ ,uu'fitu^ ^ ( 16 ) 

nil' 
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TABLE II. Ionization potentials (eV) of the porphin molecule. The 
cc-pVDZ basis set was used unless otherwise stated. The (4e, 4o) 
and (3e, 4o) active spaces were used in the CASPT2 calculations. 


Method 

Vertial IP 

Adiabatic IP 

AIP 



6.70 

6.68 

0.02 

0.77 

MP2 

8.51 

C 


1.47 

CASPT2 

6.84 

6.65 

0.18 

0.75 

CASPT2/PBE 

6.83 

6.85 

-0.02 

0.75 

Experiment 

6.9 





“ (5^) of the radical cation at the equilibrium geometry of the neutral. 
Computed using the def2-SVP basis set. 

Due to instability we could not optimize the geometry of the radical 
cation. 

Taken from Ref.ls^ 

and the Tq amplitudes of size (hocc and rtbas are 

the numbers of occupied orbitals and basis functions, re¬ 
spectively). Further optimization and distributed-memory 
parallelization of our program are warranted and will be 
performed in the future. We will also consider the level 
shift^^ other zeroth-order Hamiltonians^^iS and multistate 
extensions i 

The debugging of parts of the code in bagel, includ¬ 
ing Z-CASSCF, was facilitated by existing implementations 
in molpro3 This work has been supported by Department 
of Energy, Basic Energy Sciences (Grant No. DE-EG02- 
13ER16398) and the Air Eorce Office of Scientific Research 
Young Investigator Program (Grant No. EA9550-15-1-0031). 
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